#
import numpy as np
import os, sys
sys.path.insert(0, '/global/homes/q/qmxp55/DESI/bgstargets/py')
from io_ import get_sweep_whole, getBGSbits, flux_to_mag
from io_ import get_random, get_isdesi, get_dict, bgsmask, get_reg, get_svfields, get_svfields_fg, gaiaAEN
from cuts import getGeoCuts, bgsbut
from QA import getStats, flow, mollweide, mycmap, plot_sysdens, overdensity, hexbin
from postages_images import postages_circle
import healpy as hp
import astropy.io.fits as fits
import fitsio
import matplotlib.pyplot as plt
from matplotlib import gridspec
import matplotlib
from astropy.coordinates import SkyCoord
import astropy.units as units
%load_ext autoreload
%autoreload 2
import warnings
warnings.filterwarnings('ignore')
#
dr = 'dr9g'
survey = 'south' #is either south (DECaLS+DES) or north (BASS/MzLS)
version = '0.1.0'
filesdir = '/global/cscratch1/sd/qmxp55/bgstargets_output/'
Nranfiles = 3
reg = 'svfields_fg'
dec_resol_ns = 32.375
if (dr[:3] == 'dr9'): Nranfiles = 20 # because the randoms for dr9d have a density of 100000 = 5000*20
pathdir = os.path.abspath(os.getcwd())+'/%s-%s_%s_%s_draft_plots' %(dr, survey, reg, version)
ispathdir = os.path.isdir(pathdir)
if not ispathdir: os.mkdir(pathdir)
if reg[:8] == 'svfields': desifootprint = False
else: desifootprint = True
# for healpy
dr8pix = '/project/projectdirs/desi/target/catalogs/dr8/0.31.1/pixweight/pixweight-dr8-0.31.1.fits'
hdr = fits.getheader(dr8pix,1)
nside,nest = hdr['hpxnside'],hdr['hpxnest']
npix = hp.nside2npix(nside)
pixarea = hp.nside2pixarea(nside,degrees=True)
#load catalogues
if dr[:3] == 'dr9': N = 1
else: N = Nranfiles
cat = np.load(filesdir+dr+'/'+version+'/'+'bgstargets-'+survey+'.npy')
cat_ex = np.load(filesdir+dr+'/'+version+'/'+'extra-'+survey+'_n256.npy')
ran = np.load(filesdir+dr+'/'+dr+'_random_N'+str(N)+'.npy')
ran_ex = np.load(filesdir+dr+'/'+'extra_random_N'+str(N)+'_n256.npy')
hppix_ran = ran_ex['hppix']
ranindesi = ran_ex['desi']
#if svfields don't use desifootprint
if not desifootprint: raninreg = (ran_ex[reg]) & (ran_ex[survey])
else: raninreg = (ran_ex[reg]) & (ran_ex['desi']) & (ran_ex[survey])
hppix_cat = cat_ex['hppix']
catindesi = cat_ex['desi']
#if svfields don't use desifootprintreal amateurreal amateur
if not desifootprint: catinreg = (cat_ex[reg]) & (cat_ex[survey])
else: catinreg = catinreg = (cat_ex[reg]) & (cat_ex['desi']) & (cat_ex[survey])
#RANDOMS
keep = raninreg
ran = ran[keep]
ran_ex = ran_ex[keep]
ranindesi = ranindesi[keep]
hppix_ran = hppix_ran[keep]
raninreg = raninreg[keep]
#CATALOGUE
keepC = catinreg
cat = cat[keepC]
cat_ex = cat_ex[keepC]
catindesi = catindesi[keepC]
hppix_cat = hppix_cat[keepC]
catinreg = catinreg[keepC]
#bgsmask = bgsmask()
rancuts = getGeoCuts(ran, randoms=True)
NOTE: if using SVFIELDS footprint, use desifootprint=False otherwise the area will be less than the current SVFIELDS due desi footprint chops out part of it.
#Notes: This hpdict is used to get the area only using the randoms.
#This use the randoms within the DESI footprint and without any masking (maskrand=None) as we want the area without wholes.
hpdict0 = get_dict(cat=None, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=None,
maskrand=None, maskcat=None, Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=desifootprint, namesels=None, target_outputs=False, log=True)
#from astropy.coordinates import SkyCoord
#import astropy.units as units
#c = SkyCoord(cat['RA']*units.degree,cat['DEC']*units.degree, frame='icrs')
#galb = c.galactic.b.value # galb coordinate
#dic with default BGS selection and in DESI footprint
#Appy or not apply LG mask in randoms (rancuts['LG'])?
hpdict = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['nobs'])),
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=desifootprint, namesels={'bgs_any':20, 'bright':21, 'faint':22},
galb=cat_ex['b'], log=False, survey='bgs')
#
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
bgs_b = ((cat['BGSBITS'] & 2**(21)) != 0)
bgs_f = ((cat['BGSBITS'] & 2**(22)) != 0)
print('mean density of BGS: \t %.2f' %(hpdict['meandens_bgs_any_'+survey]))
#print('total density of BGS: \t %.2f' %(hpdict['dens_bgs_any_'+survey]))
print('Area within BGS: \t %.2f' %(hpdict['bgsarea_'+survey]))
print('No BGS galaxies in Bright and Faint: \t %i \t %i \t %i' %(
np.sum(bgs_b),
np.sum(bgs_f),
np.sum(bgs_any)
))
print('Overall density of BGS: \t %d \t %d \t %d' %(
np.round(hpdict['dens_bright_'+survey], 0),
np.round(hpdict['dens_faint_'+survey], 0),
np.round(hpdict['dens_bgs_any_'+survey], 0)
))
print('morph \t BRIGHT \t FAINT \t OVERALL')
print('------------------------------------')
for morph in list(set(cat['TYPE'][:10000])):
keep = cat['TYPE'] == morph
print('%s \t %d \t %d \t %d' %(
morph,
np.round(np.sum((keep) & (bgs_b))/hpdict['bgsarea_'+survey], 0),
np.round(np.sum((keep) & (bgs_f))/hpdict['bgsarea_'+survey], 0),
np.round(np.sum((keep) & (bgs_any))/hpdict['bgsarea_'+survey], 0),
))
if dr[:3] == 'dr9':
morphos = ['REX', 'EXP', 'DEV', 'SER']
PSF = ['PSF']
else:
morphos = ['REX ', 'EXP ', 'DEV ', 'COMP']
PSF = ['PSF ']
set(cat['TYPE'][:10000])
from QA import mollweide, mycmap, plot_sysdens
import matplotlib.pyplot as plt
from matplotlib import gridspec
import matplotlib
#
if reg[:8] == 'svfields': reg_ = reg+'_'+survey[0]
else: reg_ = reg
fig = plt.figure(figsize=(15, 15))
gs = gridspec.GridSpec(2,1)
org = 120 # centre ra for mollweide plots
projection = 'mollweide'
cm = mycmap(matplotlib.cm.jet, 10,0,1)
cmr= mycmap(matplotlib.cm.jet_r,10,0,1)
mollweide(hpdict=hpdict, namesel='bright', reg=reg_, projection=projection, n=0, org=org, cm=cm, fig=fig, gs=gs, cval=None, desifootprint=False)
mollweide(hpdict=hpdict, namesel='faint', reg=reg_, projection=projection, n=1, org=org, cm=cm, fig=fig, gs=gs, cval=None, desifootprint=False)
#mollweide(hpdict=hpdict, C=None ,namesel='any', reg='all', projection=projection, n=1, org=org, cm=cm,
# fig=fig, ws=ws, perc=(0.3,99.8), title='After linear weights', cval=(84, 2750))
file = '%s/skydens_%s_%s' %(pathdir, dr, reg)
fig.savefig(file+'.png', bbox_inches = 'tight', pad_inches = 0)
#use below code to get target densities of keep and rejected at a particular stage in flowchart...
#use hpdict0 if density computed over the total area (before applyin spatial cuts in randoms)
#use hpdict if density computed over the reduced area (after applying spatial cuts in randoms)
if False:
t = getStats(cat=cat, hpdict=hpdict0, bgsmask=bgsmask(), rancuts=rancuts, CurrentMask=['QC_FM', 'QC_FI', 'QC_FF'],
PrevMask=['BS', 'LG', 'GC', 'nobs', 'SG', 'FMC2', 'CC'],
reg=reg, regcat=catinreg, regran=raninreg)
if reg == 'south': surveylab = 'DECaLS+DES'
elif reg == 'north': surveylab = 'BASS-MzLS'
elif reg == 'desi': surveylab = 'DESI'
else: surveylab = reg
flowTitle = dr+'_'+surveylab
#this is important to read the area of the proper region in hpdict0
region = reg
order = [['BS', 'LG', 'GC'], ['nobs'], ['SG'], ['FMC2', 'CC'], ['QC_FM', 'QC_FI', 'QC_FF']]
flowNominal, _, _ = flow(cat=cat, hpdict=hpdict0, bgsmask=bgsmask(), rancuts=rancuts, order=order, reg=region,
regcat=catinreg, regran=raninreg, file='%s/flow_main_nominal_%s' %(pathdir, flowTitle), dr=flowTitle, program='main')
flowNominal
order = [ ['SG'], ['BS', 'LG', 'GC'], ['nobs'], ['QC_FM', 'QC_FI', 'QC_FF'], ['FMC2', 'CC']]
flowGalview, _, _ = flow(cat=cat, hpdict=hpdict0, bgsmask=bgsmask(), rancuts=rancuts, order=order, reg=region,
regcat=catinreg, regran=raninreg, file='%s/flow_main_galview_%s' %(pathdir,flowTitle), dr=flowTitle, program='main')
flowGalview
Firs we do the stacking around the bright stars (BS mask) using BGS but BS.
For the second stacking, we use BGS and do the stacking around the Medium Stars (MS mask).
#star catalogue
stars = np.load('/global/cscratch1/sd/qmxp55/pauline/stars_GAIA_TYCHO_13.npy')
hppix_stars = hp.ang2pix(nside,(90.-stars['DEC'])*np.pi/180.,stars['RA']*np.pi/180.,nest=True) # catalogue hp pixels array
if reg[:8] == 'svfields': starsinreg = get_svfields_fg(stars['RA'], stars['DEC'])
else: starsinreg = get_reg(reg=reg, hppix=hppix_stars)
#cstars = SkyCoord(stars['RA']*units.degree,stars['DEC']*units.degree, frame='icrs')
#galb_stars = cstars.galactic.b.value # galb coordinate
stars = stars[(starsinreg) & (stars['DEC'] > -25)]
#
#bgslist = ['LG', 'GC', 'nobs', 'SG', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF']
# get BGS objects without BS bit set on
#bgsbutBS = np.ones_like(cat['RA'], dtype=bool)
#for key in bgslist:
# keep = ((cat['BGSBITS'] & 2**(bgsmask()[key])) != 0)
# bgsbutBS &= keep
#bgsbutBS &= cat['RMAG'] < 20
bgsbutBS = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['BS'], bgsmask=bgsmask(), rlimit=20)
log = True
nbins = 150
#Dustin_radii's
mag = np.linspace(0, 20, 50)
BS_radii = []
new_BS_radii = []
for i,j in enumerate(mag):
BS_radii.append([j, np.minimum(1800., 150. * 2.5**((11. - j)/3.)) * 0.262/1])
new_BS_radii.append([j, 0.5 * np.minimum(1800., 150. * 2.5**((11. - j)/3.)) * 0.262/1])
#radius = np.minimum(1800., 150. * 2.5**((11. - mag)/3.)) * 0.262/3600.
plt.figure(figsize=(20, 20))
plt.title(r'%s %s' %(dr, survey), size=20)
plt.axis('off')
_ = overdensity(cat[bgsbutBS], stars, BS_radii, 'MAG', 35, density=False,
magbins=[10,13,4], radii_2=new_BS_radii, grid=[2,2], SR=[0., 250.], scaling=False, nbins=nbins,
SR_scaling=3,logDenRat=[-5, 5], radii_bestfit=False, annulus=None, bintype='1',
filename='%s/2d_stack_BS_%s_%s' %(pathdir, dr, survey), log=log)
plt.figure(figsize=(10,5))
plt.plot(np.transpose(BS_radii)[0], np.transpose(BS_radii)[1], lw=2, c='k', label='BS radii')
plt.plot(np.transpose(new_BS_radii)[0], np.transpose(new_BS_radii)[1], lw=2, c='r', label='NEW BS radii')
plt.legend()
plt.grid()
from QA import circular_mask_radii_func
from io_ import query_catalog_mask
BS_t = query_catalog_mask(cat['RA'], cat['DEC'], stars, BS_radii, nameMag='MAG', diff_spikes=False,
length_radii=None, widht_radii=None, return_diagnostics=True, bestfit=False, log=False)
new_BS_t = query_catalog_mask(cat['RA'], cat['DEC'], stars, new_BS_radii, nameMag='MAG', diff_spikes=False,
length_radii=None, widht_radii=None, return_diagnostics=True, bestfit=False, log=False)
#masking radius for each object cloosest stars in cat
mask_radii = circular_mask_radii_func(BS_t[1]['w1_source'], BS_radii, bestfit=False)
#masking radius reescaled for each object
mask_radii_res = BS_t[1]['d2d_source']/mask_radii
BS_diff = (BS_t[0]) & (~new_BS_t[0]) & (cat['RMAG'] < 20)
PSF = ((cat['TYPE'] == 'PSF') | (cat['TYPE'] == 'PSF '))
refcat = cat['REF_CAT']
if isinstance(np.atleast_1d(refcat)[0], str):
LG = [(rc[0] == "L") if len(rc) > 0 else False for rc in refcat]
else:
LG = [(rc.decode()[0] == "L") if len(rc) > 0 else False for rc in refcat]
print('objects within BS radii defs & r<20: \t %i (PSF) \t %i (no-PSF) \t %i (LSLGA)'
%(np.sum((BS_diff) & (PSF)), np.sum((BS_diff) & (~PSF)), np.sum((BS_diff) & (LG))))
print('objects within BS radii defs & r<20 & BGS: \t %i (PSF) \t %i (no-PSF) \t %i (LSLGA)'
%(np.sum((BS_diff) & (PSF) & (bgsbutBS)), np.sum((BS_diff) & (~PSF) & (bgsbutBS)), np.sum((BS_diff) & (LG))))
$\Delta \vec{r}/r_{BS}$
#
rows, cols = 1, 2
coord = {'r':cat['RMAG'], '$\Delta r/r_{BS}$':mask_radii_res}
area = hpdict0['bgsarea_'+survey]
#debug = cat['RMAG'] < 16
masks = {'BS_old not in BS_new & r<20':(BS_diff), 'BS_old not in BS_new & r<20 & BGS':(BS_diff) & (bgsbutBS)}
hline, vline = 0.5, None
vmin, vmax = 1, None
fig = plt.figure(figsize=(8*cols,6*rows))
gs = gridspec.GridSpec(rows, cols, hspace=0.25, wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
for i, key, val in zip(range(len(masks)), masks.keys(), masks.values()):
if (i%cols==0): ylab=True
else: ylab = False
hexbin(coord=coord, catmask=(val), n=i, bins='log', title=key, cmap=cmap,
ylab=ylab, xlab=True, vline=vline, hline=hline, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, fmcline=False,
file='%s/r_mask_radii_BSold_not_BSnew_%s' %(pathdir, survey), fracs=False, area=area, cbar='horizontal',
xlim=[12, 20.1], ylim=[0.4, 1.0])
#
veto = {'BGS':(bgs_any),'Not BGS':(~bgs_any)}
info = {'r':cat['RMAG'], 'T':cat['TYPE'].astype(str), 'ref':cat['REF_CAT'].astype(str), 'd2d':np.round(mask_radii_res, 1)}
layer = dr+'-'+survey
#'dr8-model'
#'dr8-resid'
#layer2=layer+'-model'
masks = {'rmag<19':(BS_diff) & (bgsbutBS) & (cat['RMAG'] < 19), 'rmag>19':(BS_diff) & (bgsbutBS) & (cat['RMAG'] > 19)}
for key, val in zip(masks.keys(), masks.values()):
idx = list(np.where(val))[0]
print('sample size: %i' %(len(idx)))
#if we want to compare the same objects with another catalogue (i.e., DR7, DR8) select only
#the right number of indexes as your grid to avoid the random selection in postages_circle.
postages_circle(coord=[cat['RA'], cat['DEC']], centeridx=idx, veto=veto, info=info, scale=0.262,
scale_unit='pixscale', layer=layer, radius=4/3600., m=2, grid=[12,6],
savefile='%s/%s_BSold_not_BSnew_bgs_%s' %(pathdir, layer, key), layer2=layer+'-resid', layer2Mode='merge',
isLG=False, title=None, markers=False, colorkey=True)
plt.show()
from io_ import get_msmask
#load mask sources objects from SWEEPS DR8
masksources = np.load('/global/cscratch1/sd/qmxp55/bgstargets_output/dr8/dr8_sweep_whole_maskbitsources.npy')
#get the medium bright stars
starsMS = get_msmask(masksources)
hppix_starsMS = hp.ang2pix(nside,(90.-starsMS['DEC'])*np.pi/180.,starsMS['RA']*np.pi/180.,nest=True) # catalogue hp pixels array
if reg[:8] == 'svfields': starsMSinreg = get_svfields_fg(starsMS['RA'], starsMS['DEC'])
else: starsMSinreg = get_reg(reg=reg, hppix=hppix_starsMS)
#cstars = SkyCoord(stars['RA']*units.degree,stars['DEC']*units.degree, frame='icrs')
#galb_stars = cstars.galactic.b.value # galb coordinate
starsMS = starsMS[(starsMSinreg) & (starsMS['DEC'] > -25)]
#
log = True
nbins = 150
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
plt.figure(figsize=(20, 20))
plt.title(r'%s %s' %(dr, survey), size=20)
plt.axis('off')
_ = overdensity(cat[bgs_any], starsMS, BS_radii, 'MAG', 35, density=False,
magbins=[13,16,3], radii_2=new_BS_radii, grid=[2,2], SR=[0., 110.], scaling=False, nbins=nbins,
SR_scaling=3,logDenRat=[-5, 5], radii_bestfit=False, annulus=None, bintype='2',
filename='%s/2d_stack_MS_%s_%s' %(pathdir, dr, survey), log=log)
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
## get BGS objects without BS bit set on
#bgsbutSG = np.ones_like(cat['RA'], dtype=bool)
#for key in bgslist:
# keep = ((cat['BGSBITS'] & 2**(bgsmask()[key])) != 0)
# bgsbutSG &= keep
#bgsbutSG &= cat['RMAG'] < 20
bgsbutSG = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['SG'], bgsmask=bgsmask(), rlimit=20)
finite = np.ones_like(cat['RMAG'], dtype='?')
for i in ['RMAG', 'GMAG', 'ZMAG','FLUX_R', 'RFIBERMAG']:
finite &= np.isfinite(cat[i])
Ared = hpdict['bgsarea_'+reg]
#
Grr = cat['G'] - flux_to_mag(cat['FLUX_R'])
coord = {'g-z':cat['GMAG'] - cat['ZMAG'], 'G-rr':Grr}
PSF = (cat['TYPE'] == 'PSF ') | (cat['TYPE'] == 'PSF')
#if dr[:3] == 'dr9': morphos = ['REX', 'EXP', 'DEV', 'SER']
#else: morphos = ['REX ', 'EXP ', 'DEV ', 'COMP']
rows, cols = 2, 2
#debug = cat['RMAG'] < 16
mask = (catinreg) & (bgsbutSG) & (finite) #& (debug)
hline, vline = 0.6, None
vmin, vmax = 1, None
fig = plt.figure(figsize=(8*cols,6*rows))
gs = gridspec.GridSpec(rows, cols, hspace=0.2, wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
for i, morpho in enumerate(morphos):
if (i%cols==0): ylab=True
else: ylab = False
if (morpho == PSF):
vmax, cbar = None, 'vertical'
else:
vmax, cbar = 100, 'panel'
if i < 2: xlab = False
else: xlab = True
morphomask = cat['TYPE'] == morpho
hexbin(coord=coord, catmask=((mask) & (morphomask)), n=i, bins='log', title=morpho, cmap=cmap,
ylab=ylab, xlab=xlab, vline=vline, hline=hline, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1,
file='%s/gz_Grr_bgsbutSG_hexbin_extended_%s' %(pathdir, survey), fracs=True, area=Ared, cbar=cbar)
fig = plt.figure(figsize=(8*1,8*1))
gs = gridspec.GridSpec(1, 1, hspace=0.15, wspace=0.10)
hexbin(coord=coord, catmask=((mask) & (PSF)), n=0, bins='log', title='PSF', cmap=cmap,
ylab=True, vline=vline, hline=hline, fig=fig, gs=gs, vmin=None, vmax=1000, mincnt=1,
file='%s/gz_Grr_bgsbutSG_hexbin_psf_%s' %(pathdir, survey), fracs=True, area=Ared, cbar='horizontal')
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'SG', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
## get BGS objects without BS bit set on
#bgsbutFMC2 = np.ones_like(cat['RA'], dtype=bool)
#for key in bgslist:
# keep = ((cat['BGSBITS'] & 2**(bgsmask()[key])) != 0)
# bgsbutFMC2 &= keep
#bgsbutFMC2 &= cat['RMAG'] < 20
bgsbutFMC2 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['FMC2'], bgsmask=bgsmask(), rlimit=20)
#
rows, cols = 2, 2
coord = {'r':cat['RMAG'], 'rfibmag':cat['RFIBERMAG']}
#debug = cat['RMAG'] < 16
mask = (catinreg) & (bgsbutFMC2) & (finite) #& (debug)
hline, vline = None, None
vmin, vmax = 1, 200
fig = plt.figure(figsize=(8*cols,6*rows))
gs = gridspec.GridSpec(rows, cols, hspace=0.25, wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
for i, morpho in enumerate(morphos):
if (i%cols==0): ylab=True
else: ylab = False
if morpho == PSF:
vmax, cbar = None, 'vertical'
else:
vmax, cbar = 2000, 'panel'
if i < 2: xlab = False
else: xlab = True
morphomask = cat['TYPE'] == morpho
hexbin(coord=coord, catmask=((mask) & (morphomask)), n=i, bins='log', title=morpho, cmap=cmap,
ylab=ylab, xlab=xlab, vline=vline, hline=hline, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, fmcline=True,
file='%s/r_rfibmag_bgsbutFMC2_hexbin_extended_%s' %(pathdir, survey), fracs=False, area=Ared, cbar=cbar)
fig = plt.figure(figsize=(8*1,8*1))
gs = gridspec.GridSpec(1, 1, hspace=0.15, wspace=0.10)
hexbin(coord=coord, catmask=((mask) & (PSF)), n=0, bins='log', title='PSF', cmap=cmap,
ylab=True, vline=vline, hline=hline, fig=fig, gs=gs, vmin=None, vmax=1000, mincnt=1, fmcline=True,
file='%s/r_rfibmag_bgsbutFMC2_hexbin_psf_%s' %(pathdir, survey), fracs=False, area=Ared, cbar='horizontal')
bgsbutCC = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['CC'], bgsmask=bgsmask(), rlimit=20)
#
import matplotlib.patches as patches
coord = {'g-r':cat['GMAG'] - cat['RMAG'], 'r-z':cat['RMAG'] - cat['ZMAG']}
mask = (catinreg) & (bgsbutCC) & (finite) # & (cat['RMAG'] < 14)
fig = plt.figure(figsize=(10,10))
gs = gridspec.GridSpec(1, 1, hspace=0.15, wspace=0.10)
hexbin(coord=coord, catmask=mask, n=0, bins='log', title=None, cmap=cmap,
ylab=True, vline=None, hline=None, fig=fig, gs=gs, xlim=(-2, 7), ylim=(-4, 5),
vmin=None, vmax=None, mincnt=1, fmcline=False,
file=None, fracs=False, area=Ared, cbar='horizontal')
plt.plot(np.linspace(-1, 4, 3), np.full(3, -1), lw=2, c='r')
plt.plot(np.linspace(-1, 4, 3), np.full(3, 4), lw=2, c='r')
plt.plot(np.full(3, -1), np.linspace(-1, 4, 3), lw=2, c='r')
plt.plot(np.full(3, 4), np.linspace(-1, 4, 3), lw=2, c='r')
file='%s/gr_rz_bgsbutCC_hexbin_%s' %(pathdir, survey)
fig.savefig(file, bbox_inches = 'tight', pad_inches = 0)
from QA import plot_venn3
# get BGS but QCs
bgsbutQCs = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM', 'QC_FI', 'QC_FF'], bgsmask=bgsmask(), rlimit=20)
FM = ((cat['BGSBITS'] & 2**(11)) == 0) #rejects by FM
FI = ((cat['BGSBITS'] & 2**(12)) == 0) #rejects by FI
FF = ((cat['BGSBITS'] & 2**(13)) == 0) #rejects by FF
area = hpdict0['bgsarea_'+survey]
#
keepB = (bgsbutQCs) & (cat['RMAG'] < 19.5)
keepF = (bgsbutQCs) & (cat['RMAG'] > 19.5) & (cat['RMAG'] < 20)
plot_venn3(A=(keepB & FM), B=(keepB & FI), C=(keepB & FF), norm=area,
labels=['FM', 'FI', 'FF'], file=pathdir+'/venn_QCs_bright_%s' %(survey), title='BRIGHT', colors = ['deepskyblue', 'springgreen', 'salmon'])
plot_venn3(A=(keepF & FF), B=(keepF & FI), C=(keepF & FM), norm=area,
labels=['FF', 'FI', 'FM'], file=pathdir+'/venn_QCs_faint_%s' %(survey), title='FAINT', colors = ['salmon' ,'springgreen', 'deepskyblue'])
Tractor force to fit as PSF the GAIA objects that satisfy the following conditions:
For $G<19$: astrometric_excess_noise $< 10^{0.5}$,
For $G≥19$: astrometric_excess_noise $ < 10^{0.5+0.2(G−19)}$.
G is the GAIA photometric $G$-band and the astrometric excess noise ($AEN$) is the... Those PSF that are in GAIA are treated as Stars according to TRACTOR and the GAIA objects that do not fulfill that condition are treated as GALAXIES. This is different to the way we do star-galaxy separation for BGS, we use GAIA but rather than look at the $AEN$ we look at the colour distribution $G-rr$.
The PSF GAIA objects (from above definition) are treated differently than other PSF objects and they have a local fitting, the problem is that maybe those are not stars and might be galaxies instead that photometry is being altered by this local fitting.
Here we are going to find out what are the differences between our star-galaxy classification and the one above for GAIA objects. We're going to:
TRACTOR PSF and TRACTOR non-PSF:inGAIA = (cat['REF_CAT'] == 'G2') & (cat['RMAG'] < 20)
grr_gal = ((cat['BGSBITS'] & 2**(6)) != 0) & (cat['RMAG'] < 20)
grr_stars = ((cat['BGSBITS'] & 2**(6)) == 0) & (cat['RMAG'] < 20)
AEN_star, AEN_gal = gaiaAEN(inGAIA=inGAIA, size=len(cat['RA']), G=cat['G'], AEN=cat['AEN'], dr=dr)
PSF = (cat['TYPE'] == 'PSF ') | (cat['TYPE'] == 'PSF')
#inGAIA2 = (cat['G'] != 0)
#
coord = {'G':cat['G'], r'$\log(AEN)$':np.log10(cat['AEN'])}
#masks = [(grr_stars) & (inGAIA) & (cat['RMAG'] < 14), (grr_gal) & (inGAIA) & (cat['RMAG'] < 14)]
masks = [(grr_stars) & (inGAIA), (grr_gal) & (inGAIA)]
titles = ['stars with G-rr (r < 20)', 'galaxies with G-rr (r < 20)']
vmin, vmax = 1, None
fig = plt.figure(figsize=(9*len(masks), 7))
gs = gridspec.GridSpec(1, len(masks),hspace=0.1,wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
xlim=[10, 22]
ylim=[-2, 1.5]
for i, mask in enumerate(masks):
#if (i%len(masks)==0): xlab=True
#else: xlab = False
#if morpho == 'PSF ': vmax = None
#morphomask = cat['TYPE'] == morpho
hexbin(coord=coord, catmask=mask, n=i, bins='log', title=titles[i], cmap=cmap, xlab=True,
ylab=True, vline=None, hline=None, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, fmcline=False,
file=None, area=None, xlim=xlim, ylim=ylim, cbar='horizontal')
#if i == 0:
# x, y = coord.keys()
# mask1 = (gaia_star) & (dr8_gama) & (dr8_z > 0.002) & (bgs)
# mask2 = (gaia_star) & (dr8_gama) & (dr8_z > 0.002) & (~bgs)
# plt.scatter(coord[x][mask2], coord[y][mask2], s=2, c='royalblue', label=r'Not in BGS')
# plt.scatter(coord[x][mask1], coord[y][mask1], s=2, c='k', label=r'in BGS')
# plt.legend()
x_ = np.linspace(19, xlim[1], 4)
plt.plot(np.linspace(xlim[0], 18, 3), np.full(3, 0.5), color='r', ls='-', lw=2)
plt.plot(np.full(3, 18), np.linspace(ylim[0], 0.5, 3), color='r', ls='-', lw=2)
plt.plot(x_, 0.5 + 0.2*(x_ - 19.), color='r', ls='--', lw=2, alpha=0.7)
plt.plot(np.linspace(18,19,3), np.full(3, 0.5), color='r', ls='--', lw=2, alpha=0.7)
file = pathdir+'/G_logAEN'
fig.savefig(file+'.png', bbox_inches = 'tight', pad_inches = 0)
#
area = hpdict0['bgsarea_'+survey]
BS = ((cat['BGSBITS'] & 2**(0)) != 0)
gal_in_aenstar = (inGAIA) & (grr_gal) & (cat['G'] < 18)*(np.log10(cat['AEN']) < 0.5)
print('G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5): \t %i, \t %.3f (1/deg^2)'
%(np.sum(gal_in_aenstar), np.sum(gal_in_aenstar)/area))
print('G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5) & BS: \t %i, \t %.3f (1/deg^2)'
%(np.sum(gal_in_aenstar & (BS)), np.sum(gal_in_aenstar & (BS))/area))
print('G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5) & BGS: \t %i, \t %.3f (1/deg^2)'
%(np.sum(gal_in_aenstar & (bgs_any)), np.sum(gal_in_aenstar & (bgs_any))/area))
#
veto = {'BGS':(bgs_any),'Not BGS':(~bgs_any)}
info = {'r':cat['RMAG'], 'T':cat['TYPE'].astype(str), 'ref':cat['REF_CAT'].astype(str)}
layer = dr+'-'+survey
#'dr8-model'
#'dr8-resid'
#layer2=layer+'-model'
idx = list(np.where((gal_in_aenstar) & (bgs_any)))[0]
print('sample size: %i' %(len(idx)))
#if we want to compare the same objects with another catalogue (i.e., DR7, DR8) select only
#the right number of indexes as your grid to avoid the random selection in postages_circle.
postages_circle(coord=[cat['RA'], cat['DEC']], centeridx=idx, veto=veto, info=info, scale=0.262,
scale_unit='pixscale', layer=layer, radius=4/3600., m=2, grid=[4,5],
savefile='%s/%s_Grr_gal_AEN_stars' %(pathdir, layer), layer2=layer+'-resid', layer2Mode='merge',
isLG=False, title=None, markers=False, colorkey=True)
#
fig = plt.figure(figsize=(18,12))
log = True
Grr = cat['G'] - flux_to_mag(cat['FLUX_R'])
classes = {'stars_Grr':(grr_stars) & (inGAIA), 'gal_Grr':(grr_gal) & (inGAIA), 'stars_aen': (AEN_star) & (cat['RMAG'] < 20), 'gal_aen':(AEN_gal) & (cat['RMAG'] < 20),
'stars_Grr_PSF':(grr_stars) & (inGAIA) & (PSF), 'gal_Grr_PSF':(grr_gal) & (inGAIA) & (PSF),
'stars_aen_PSF':(AEN_star) & (PSF) & (cat['RMAG'] < 20), 'gal_aen_PSF':(AEN_gal) & (PSF) & (cat['RMAG'] < 20),
'stars_Grr_noPSF':(grr_stars) & (inGAIA) & (~PSF), 'gal_Grr_noPSF':(grr_gal) & (inGAIA) & (~PSF),
'stars_aen_noPSF':(AEN_star) & (~PSF) & (cat['RMAG'] < 20), 'gal_aen_noPSF':(AEN_gal) & (~PSF) & (cat['RMAG'] < 20)}
for i, j in enumerate(['Grr', 'aen', 'Grr_PSF', 'aen_PSF', 'Grr_noPSF', 'aen_noPSF']):
XS = (classes['stars_' + j])
XG = (classes['gal_' + j])
tot_gaia = np.sum((inGAIA))
plt.subplot(3, 2, i+1)
if i == 0:
plt.title(r'$G-rr$ classification (r < 20)', size=22)
if i == 1:
plt.title(r'$AEN$ classification (r < 20)', size=22)
bins = np.linspace(-2.1, 5, 60)
plt.hist(Grr[XS], bins=bins, log=log, alpha=0.7, label='%s stars, $f$=%2.3g' %(j, np.sum(XS)/tot_gaia))
plt.hist(Grr[XG], bins=bins, log=log, alpha=0.4, label='%s galaxies, $f$=%2.3g' %(j, np.sum(XG)/tot_gaia))
if i == 0:
plt.axvline(0.6, ls='--', lw=2, c='k', label='SG threshold, Grr = 0.6')
plt.text(-2.1, 2*10**7, '$N_{tot}$=%s' %(tot_gaia))
plt.axvline(0.6, ls='--', lw=2, c='k')
if (i == 4) or (i == 5):
plt.xlabel(r'$G-rr$', size=20)
if i == 2:
plt.ylabel(r'$N$', size=20)
plt.legend()
plt.show()
filename = '%s/gaia_Grr_comparison_with_aen' %(pathdir)
fig.savefig(filename+'.png') #bbox_inches = 'tight', pad_inches = 0
from matplotlib_venn import venn2, venn2_circles
fig = plt.figure(figsize=(5,5))
sf = 2
area = hpdict0['bgsarea_'+survey]
a = classes['gal_Grr']
b = classes['gal_aen']
c = (a) & (b)
a1 = round(((np.sum(a) - np.sum(c))/area), sf)
b1 = round(((np.sum(b) - np.sum(c))/area), sf)
c1 = round(np.sum(c)/area, sf)
plt.title(r'GALAXIES', size=18)
labels = ('Grr', 'AEN')
venn2([a1, b1, c1], set_labels = labels)
c=venn2_circles([a1, b1, c1], linestyle='solid', linewidth=1, color="k")
filename = '%s/venn_gaia_Grr_comparison_with_aen' %(pathdir)
fig.savefig(filename+'.png', bbox_inches = 'tight', pad_inches = 0)
#
if reg[:8] == 'svfields': reg_ = reg+'_'+survey[0]
else: reg_ = reg
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi = (hpdict['isdesi']) & (hpdict['bgsfracarea']>0) #this should be rename as "mainreg"
nx = 20
b0, m0 = plot_sysdens(hpdicttmp=hpdict, namesels=['bgs_any'], regs=[reg_], syst='stardens', mainreg=isdesi,
xlim=None, n=0, nx=nx, clip=True, denslims=False, ylab=True, weights=True,
fig=fig, gs=gs, title='test', label=True)
#dic WEIGHTED
ws = 1./((m0)*hpdict['stardens'] + b0)
hpdict_ws = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels={'bgs_any':20, 'bright':21, 'faint':22},
galb=cat_ex['b'], log=True, survey='bgs', ws=ws)
from QA import pixhistregs
# dr8_south+north : density distributions + systematics
# settings
isdesi = (hpdict['isdesi']) & (hpdict['bgsfracarea']>0) #this should be rename as "mainreg"
systs = ['log10_stardens','ebv','psfsize_g', 'psfsize_r', 'psfsize_z','galdepth_g', 'galdepth_r', 'galdepth_z']
nx = 20
cols = ['0.5','b','g','r']
#regs = [reg]
namesels=['bgs_any', 'bright', 'faint']
#
# looping on subselections
#for key, title in zip([None, ws], ['UNWEIGHTED', 'WEIGHTED']):
for key, title in zip([hpdict, hpdict_ws], ['UNWEIGHTED', 'WEIGHTED']):
# systematics
fig = plt.figure(figsize=(20,15))
gs = gridspec.GridSpec(3,3,hspace=0.30,wspace=0.25)
for i in range(len(systs)+1):
if i == 0:
axinfo = fig.add_subplot(gs[i])
# infos
axinfo.axis('off')
#axinfo.text(0.5,0.9,namesel,fontsize=20,fontweight='bold',ha='center',transform=axinfo.transAxes)
tmpy = 0.7
#for regi,col in zip(regs,cols):
#if reg == 'south':reg_ = 'DECaLS+DES'
#else: reg_ = reg
for namesel,col in zip(namesels,cols):
tmpstr = surveylab+' & '+namesel+' : '+r'$\overline{n}=$'+'%.0f'%key['meandens_'+namesel+'_'+reg]+r' deg$^{-2}$'
axinfo.text(0.05,tmpy,tmpstr,color=col,fontsize=15,transform=axinfo.transAxes)
tmpy -= 0.15
else:
syst = systs[i-1]
if (i%3==0) or (i==1): ylab=True
else: ylab = False
if i == 1: label = True
else: label = False
#if title == 'UNWEIGHTED':
# weights = False
# onlyweights=False
#else:
# weights = True
# onlyweights=True
plot_sysdens(hpdicttmp=key, namesels=namesels, regs=[reg], syst=syst, mainreg=isdesi, xlim=None, n=i, nx=nx, clip=True,
denslims=False, ylab=ylab, weights=False, fig=fig, gs=gs, label=label, title=title,
ws=None)
# save fig
fig.savefig('%s/systematics_main_bgs_%s_%s.png' %(pathdir, survey, title), bbox_inches = 'tight', pad_inches = 0)
print('')
print('')
namesels = {}
for i in ['COMP', 'DEV ', 'EXP ', 'PSF ', 'REX ']:
namesels[i] = (cat['TYPE'] == i) & (bgs_any)
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels=namesels, galb=cat_ex['b'], log=True, survey='custom')
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['COMP', 'DEV ', 'EXP ', 'REX '], regs=[reg],
syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test2', label=True)
namesels = {}
mags = np.linspace(16, 20, 9)
for i in range(len(mags[:-1])):
namesels['%s_%s' %(str(mags[i]), str(mags[i+1]))] = (cat['RMAG'] > mags[i]) & (cat['RMAG'] < mags[i+1]) & (bgs_any)
namesels.keys()
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels=namesels, galb=cat_ex['b'], log=True, survey='custom')
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['16.0_16.5', '16.5_17.0', '17.0_17.5', '17.5_18.0'], regs=[reg],
syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test3', label=True)
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['18.0_18.5', '18.5_19.0', '19.0_19.5', '19.5_20.0'], regs=[reg],
syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test4', label=True)
Below code only needs the randoms.
#
from QA import set_mwd, get_radec_mw
from io_ import get_isdes
# dict without geometrical mask and in WHOLE SWEEPS footprint
hpdict0 = get_dict(cat=None, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=None,
maskrand=None, maskcat=None, Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=False, namesels=None, target_outputs=False, log=False)
# dict without geometrical mask and in DESI footprint
hpdict1 = get_dict(cat=None, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=None,
maskrand=None, maskcat=None, Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels=None, target_outputs=False, log=False)
#
for i in ['area_all', 'bgsarea_south', 'bgsarea_north', 'bgsarea_des']:
print('%s: \t %.3f(whole) \t %.3f(desifootprint)' %(i, hpdict0[i], hpdict1[i]))
if survey == 'south':
# DECaLS area in the NGC/SGC within DESI footprint
area_decals_ngc = hpdict1['bgsfracarea'][(hpdict1['issouth']) & (hpdict1['galb'] > 0)].sum() * pixarea
area_decals_sgc = hpdict1['bgsfracarea'][(hpdict1['issouth']) & (hpdict1['galb'] < 0)].sum() * pixarea
print('NGC DECaLS in DESI footprint: \t %.3f' %(area_decals_ngc))
print('SGC DECaLS in DESI footprint: \t %.3f' %(area_decals_sgc))
if survey == 'north':
# DECaLS area in the NGC/SGC within DESI footprint
area_bass_ngc = hpdict1['bgsfracarea'][(hpdict1['isnorth']) & (hpdict1['galb'] > 0)].sum() * pixarea
print('NGC BASS/MzLS in DESI footprint: \t %.3f' %(area_bass_ngc))
desitiles = fitsio.read('/global/homes/q/qmxp55/desi-tiles-viewer.fits')
org = 120 # centre ra for mollweide plots
projection = 'mollweide'
# plotting bgsfracarea
fig = plt.figure(figsize=(14, 7))
gs = gridspec.GridSpec(1,1,wspace=0.15,hspace=0)
ax = plt.subplot(gs[0],projection=projection)
_ = set_mwd(ax,org=org)
ramw,decmw = get_radec_mw(hpdict0['ra'],hpdict0['dec'],org)
tmp = (hpdict0['bgsfracarea']>0) & (hpdict0['is'+survey])
if survey == 'south':
SC = ax.scatter(ramw[tmp],decmw[tmp],s=0.4, color='gray', label=r'DECam imaging')
if survey == 'north':
SC = ax.scatter(ramw[tmp],decmw[tmp],s=0.4, color='gray', label=r'BASS/MzLS imaging')
ramw,decmw = get_radec_mw(desitiles['ra'],desitiles['dec'],org)
if survey == 'south':
tmp = (desitiles['in_desi'] == 1) & (desitiles['dec'] < dec_resol_ns)
SC = ax.scatter(ramw[tmp],decmw[tmp], marker='.', s = 220, facecolors='none', edgecolors='blue', label=r'DECam LS dr8 in DESI')
if survey == 'north':
tmp = (desitiles['in_desi'] == 1) & (desitiles['dec'] > dec_resol_ns)
SC = ax.scatter(ramw[tmp],decmw[tmp], marker='.', s = 220, facecolors='none', edgecolors='blue', label=r'BASS/MzLS dr8 in DESI')
if survey == 'south':
desitilesindes = get_isdes(desitiles['ra'],desitiles['dec'])
tmp = (desitiles['in_desi'] == 1) & (desitiles['dec'] < dec_resol_ns) & (desitilesindes)
SC = ax.scatter(ramw[tmp],decmw[tmp], marker='.', s = 220, facecolors='none', edgecolors='red', label=r'DECam DES in DESI')
gc = SkyCoord(l=np.linspace(0, 360, 50)*units.degree, b=np.full(50, 0)*units.degree, frame='galactic')
ramw,decmw = get_radec_mw(gc.icrs.ra.value, gc.icrs.dec.value,org)
SC = ax.scatter(ramw, decmw, color='r', marker='.', lw=2)
lgnd = ax.legend()
[handle.set_sizes([120.0]) for handle in lgnd.legendHandles]
fig.savefig('%s/skyplot_%s.png' %(pathdir, survey), bbox_inches = 'tight', pad_inches = 0)
# get the area covered by the number of passes from 1 to 3 and per band and for the joint bands.
Apasses = {}
for i in [0, 1, 2]:
for j in ['G', 'R', 'Z', 'all']:
if j == 'all': mask = (ran['NOBS_G'] > i) & (ran['NOBS_R'] > i) & (ran['NOBS_Z'] > i)
else: mask = (ran['NOBS_'+j] > i)
hpdicttest = get_dict(pixmapfile=dr8pix, hppix_ran=hppix_ran, maskrand=mask, Nranfiles=Nranfiles,
ranindesi=ranindesi, desifootprint=True, target_outputs=False)
for k in ['south', 'north', 'decals', 'des']:
Apasses[j+str(i+1)+'_'+k] = hpdicttest['bgsarea_'+k]
Apasses[j+str(i+1)+'_'+'all'] = hpdicttest['area_all']
print('%s >= %i: \t south: %2.3g, north: %2.3g, decals: %2.3g, des: %2.3g, all: %2.3g'
%(j, i+1, hpdicttest['bgsarea_south'], hpdicttest['bgsarea_north'], hpdicttest['bgsarea_decals'],
hpdicttest['bgsarea_des'], hpdicttest['area_all']))
#
print('BAND \t >= 1 \t\t >= 2 \t\t >= 3')
print('------------------------------------------')
for band in ['G', 'R', 'Z', 'all']:
lab1 = band+'1'+'_south'
lab2 = band+'2'+'_south'
lab3 = band+'3'+'_south'
print('%s: \t %.0f \t %.0f \t %.0f' %(band, Apasses[lab1], Apasses[lab2], Apasses[lab3]))
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
bgsbut20 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20)
bgsbut201 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20.1)
bgsbut205 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20.5)
bgsbut207 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20.7)
#bgsbutCC = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['CC'], bgsmask=bgsmask())
#bgsbutFMC2 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['FMC2'], bgsmask=bgsmask())
#bgsbutBSLG = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['BS', 'LG'], bgsmask=bgsmask())
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
namesels = {}
#namesels['bgs_any'] = bgs_any
namesels['bgsbut20'] = bgsbut20
namesels['bgsbut201'] = bgsbut201
namesels['bgsbut205'] = bgsbut205
namesels['bgsbut207'] = bgsbut207
#dic with default BGS selection and in DESI footprint
#maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs']))
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=None,
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels=namesels, galb=galb, log=True, survey='custom')
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgsbut20', 'bgsbut201', 'bgsbut205', 'bgsbut207'], regs=[reg],
syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
#plot_sysdens(hpdicttmp=hpdict_tmp, namesels=[ 'bgs_19.0_19.5', 'bgs_19.5_20.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=1, nx=nx, clip=True,
# denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
bgsbut20 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20)
bgsbut201 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20.1)
bgsbut205 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20.5)
bgsbut207 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20.7)
namesels = {}
namesels['bgsbut20'] = bgsbut20
namesels['bgsbut201'] = bgsbut201
namesels['bgsbut205'] = bgsbut205
namesels['bgsbut207'] = bgsbut207
#dic with default BGS selection and in DESI footprint
#maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs']))
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels=namesels, galb=galb, log=True, survey='custom')
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgsbut20', 'bgsbut201', 'bgsbut205', 'bgsbut207'], regs=[reg],
syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test2', label=True)
#plot_sysdens(hpdicttmp=hpdict_tmp, namesels=[ 'bgs_19.0_19.5', 'bgs_19.5_20.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=1, nx=nx, clip=True,
# denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
#
ws = 1./((-3.96*10**(-5))*hpdict['stardens'] + 1.03)
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
plot_sysdens(hpdicttmp=hpdict, namesels=['bgs_any', 'bright', 'faint'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=nx, clip=True,
denslims=False, ylab=True, weights=True, fig=fig, gs=gs, title='test', label=True, ws=ws, onlyweights=True)
hpdict['meandens_bgs_any_south']
hpdict_tmp['meandens_bgs_any_south']
#
#ws = 1./((-3.96*10**(-5))*hpdict['stardens'] + 1.03)
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgs_any', 'bright', 'faint'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=nx, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True, ws=None, onlyweights=False)
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
namesels = {}
mags = [15, 18, 18.5, 19.0, 19.5, 20.0]
for num,i in enumerate(mags[:-1]):
namesels['bgs_%s_%s' %(str(mags[num]), str(mags[num+1]))] = (bgs_any) & (cat['RMAG'] > mags[num]) & (cat['RMAG'] < mags[num+1])
#dic with default BGS selection and in DESI footprint
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat,
maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
maskcat=None,
Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None,
desifootprint=True, namesels=namesels, galb=galb, log=True, survey='custom')
fig = plt.figure(figsize=(18,5))
gs = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgs_15_18', 'bgs_18_18.5', 'bgs_18.5_19.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=nx, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
plot_sysdens(hpdicttmp=hpdict_tmp, namesels=[ 'bgs_19.0_19.5', 'bgs_19.5_20.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=1, nx=nx, clip=True,
denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
bgs_bright = ((cat['BGSBITS'] & 2**(21)) != 0)
bgs_faint = ((cat['BGSBITS'] & 2**(22)) != 0)
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)